#!/bin/csh
##############################################################################
#job file to produce postscript graphic files of TISC results with GMT 4.0
#	Daniel Garcia-Castellanos, 2007.
##############################################################################
# -This script uses bc unix command as a calculator. 
##############################################################################

set prj 	= $1

set width = 8

source $tisc_dir/script/tisc.common.gmt.job

#TOPOGRAPHY+DRAINAGE:
#grdgradient $tmp.top.grd.tmp -A0 -G$tmp.intensity.top.grd.tmp  -Nt1
psbasemap -JX$size -R$Region -B$ticklabels\:"x, longitude (km)":/$ticklabels\:"y, latitude (km)":NseW \
	-Y10 -X2.5 -K -P > $ps

source $tisc_dir/script/tisc.common.topo+drainage.gmt.job

#PRECIPITATION
if (-r $prj.xyw) then 
	awk '{print($1,$2,$10/1e3)}' $prj.xyw |\
		xyz2grd	-G$tmp.rain.grd.tmp -I$dx/$dy -R -H2
	grdclip $tmp.rain.grd.tmp -G$tmp.rain.grd.tmp -Sb$min_rain/NaN
	grdimage $tmp.rain.grd.tmp -C$tmp.rain_cpt.tmp -JX -R -Q -O -K >> $ps
	grdcontour $tmp.top.grd.tmp -JX -R -C$tmp.topo_palet.tmp -W1 -O -K >> $ps
endif

if (-r $prj.ice) then 
#ICE
#awk '{if ($11>0) print $1, $2}' $prj.xyw | \
#	psxy -JX -R -Ss$size_ice_squares -G255 -O -K >> $ps 
awk '{if ($4>2) print($1,$2,$4)}' $prj.ice |\
	xyz2grd	-G$tmp.ice.grd.tmp -I$dx/$dy -R -H2
awk '{if ($4>0) print($1,$2,$3+$4)}' $prj.ice |\
	xyz2grd	-G$tmp.icetop.grd.tmp -I$dx/$dy -R -H2
grdgradient $tmp.icetop.grd.tmp -A0 -G$tmp.intensity.icetop.grd.tmp  #-Nt1
#.0004 MUL is normal for large scale
grdmath $tmp.intensity.icetop.grd.tmp .003 MUL = $tmp.intensity.icetop.grd.tmp
grdview $tmp.ice.grd.tmp -C$tmp.ice_cpt.tmp  -JX -R \
	-I$tmp.intensity.icetop.grd.tmp -Ts -K -O >> $ps
awk '{print $1, $2, $4}' $prj.ice | \
	xyz2grd -G$tmp.ice_thick.grd.tmp -I$dx/$dy -R -H2
grdcontour $tmp.ice_thick.grd.tmp -JX -R -C100 -A500 -Wa4/50/50/255 -Wc1/50/50/255 -O -K >> $ps
grdcontour $tmp.icetop.grd.tmp -JX -R -C$tmp.topo_palet.tmp -W1 -A -O -K >> $ps
awk '{if (NR>2 && $4>2) {a=$5+$7;b=$6+$8;l=sqrt(a*a+b*b); if (l>00) print $1,$2, atan2(b,a)*180/3.1415927, l/500}}'  $prj.ice | \
	psxy -JX -R -Svb.01/.4/.2n2 -W1 -G0 -O -K >> $ps 
endif

pstext	-JX -R -O -G0 -K <<END >> $ps 
	$xtitle	$ytitle	12 0 1 3	topography & drainage
END
pstext	-JX -R -G0 -W255/255/255 -O -K <<END >> $ps 
	$xtime	$ytime	13 0 1 9	$Timenow Myr
END
#psscale -C$tmp.topo_palet.tmp -D6/-5/12/.3h \
#	-B/:"elevation (m)": -L -O -K >> $ps	

if (-r $prj.INT_ROSE) psxy $prj.INT_ROSE -JX -R -L -W9/255/0/0ta -K -O >> $ps
if (-r $prj.pfl) psxy $prj.pfl -JX -R -M -W$col_lin -H2 -O -K >> $ps 
if (-r $prj.CMP) psxy $prj.CMP -JX -R -M -W4/0 -O -K >> $ps 
#if (-r $prj.CMP.grd) grdcontour $prj.CMP.grd -JX -R -C$tmp.topo_palet.tmp -W1/200/0/0 -O -K >> $ps 
if (-r $prj.CMP.grd) grdcontour $prj.CMP.grd -JX -R -C2000 -W3/200/0/0 -L0/6000 -O -K >> $ps 
if (-r $prj.CMP.grd) grdcontour $prj.CMP.grd -JX -R -C10000 -W5/200/0/0 -O -K >> $ps 



#SECOND PANEL:
psbasemap -JX$size -R -B$ticklabels\:"x, longitude (km)":/$ticklabels\:"y, latitude (km)":Nsew \
	-X$horz_shift -K -O >> $ps
if (-r $prj.st) then 
#EROSION MINUS SEDIMENTATION RATE :
grdimage $tmp.erosrate.grd.tmp -C$tmp.erosrate_cpt.tmp -JX -R -Q -O -K >> $ps
#SEDIMENTATION MINUS EROSION:
#grdimage $tmp.st.grd.tmp -C$tmp.eros_cumul_cpt.tmp -JX -R -Q -O -K >> $ps
#grdcontour $tmp.st.grd.tmp -JX -R -C2 -L-14/-2 -W1 -O -K >> $ps
#grdcontour $tmp.st.grd.tmp -JX -R -C2 -L2/14 -W1 -O -K >> $ps
grdcontour $tmp.top.grd.tmp -JX -R -C10000 -W3/0 -O -K >> $ps
else
#GEOLOGY/LITHOLOGY:
awk '{if (NR==2) for (i=3; i<=NF; i++) dens[i]=$i; if (NR>2) {for (i=NF; i>3; i--) if ($i-$(i-1)>10) break; if ($NF>=0) print $1,$2,dens[i]; else print $1,$2,1000}}' $prj.hrz |\
	xyz2grd	-G$tmp.geol.grd.tmp -I$dx/$dy -R
grdclip $tmp.geol.grd.tmp -G$tmp.geol.grd.tmp -Sb1200/NaN
grdimage $tmp.geol.grd.tmp -C$tmp.geol_palet.tmp -Q -JX -R \
	#-I$tmp.intensity.top.grd.tmp \
	-K -O >> $ps
grdcontour $tmp.top.grd.tmp -JX -R -C10000 -W2/0 -O -K >> $ps
endif
grdcontour $tmp.top.grd.tmp -JX -R -C$tmp.topo_palet.tmp -W1 -O -K >> $ps
#if (-r $prj.rain.grd.tmp) grdcontour $tmp.rain.grd.tmp -JX -R -C.200 -W1/160/160/255 -O -K >> $ps #-T:HL
if (-r $prj.rain.grd.tmp) grdcontour $tmp.rain.grd.tmp -JX -R -C1 -W1/220/220/255 -O -K >> $ps #-T:HL

#DRAINAGE (proportional to sed. load):
if (-r $prj.xyw) then 
if (-r $prj.RIV) awk '{if (substr($0,1,1)!=">" && substr($0,1,1)!="#") print $1/1000,$2/1000; else print $0}' \
	$prj.RIV | psxy -JX -R -M -W4/70/0/0 -O -K >> $ps 
awk -v sea_level=$sea_level '{if ($4>3 && NR>2 && $5!="S" && $5!="L" && $6>=sea_level) {print $1, $2; print $7, $8; print ">"}}' $prj.xyw | \
	psxy -JX -R -M -W1/$col_masstr -H2 -O -K >> $ps 

awk -v sea_level=$sea_level '{if ($4>30 && NR>2 && $5!="S" && $5!="L" && $6>=sea_level) {print $1, $2; print $7, $8; print ">"}}' $prj.xyw | \
	psxy -JX -R -M -W3/$col_masstr -H2 -O -K >> $ps 

awk -v sea_level=$sea_level '{if ($4>300 && NR>2 && $5!="S" && $5!="L" && $6>=sea_level) {print $1, $2; print $7, $8; print ">"}}' $prj.xyw | \
	psxy -JX -R -M -W8/$col_masstr -H2 -O -K >> $ps 
endif

if (-r $prj.pfl) psxy $prj.pfl -JX -R -M -W$col_lin -H2 -O -K >> $ps 
if (-r $prj.CMP) psxy $prj.CMP -JX -R -M -W4/0/0/0 -O -K >> $ps
#if (-r $prj.INT) psxy $prj.INT -JX -R -M -W1/0/255/0 -L -O -K >> $ps
#if (-r $prj.INT2) psxy $prj.INT2 -JX -R -M -W1/0/255/0 -L -O -K >> $ps
#if (-r $prj.INT3) psxy $prj.INT3 -JX -R -M -W1/0/255/0 -L -O -K >> $ps
#if (-r $prj.INT4) psxy $prj.INT4 -JX -R -M -W1/0/255/0 -L -O -K >> $ps
pstext	-JX -R -O -G0 -K <<END >> $ps
#	$xtitle	$ytitle	12 0 1 3	cumulative sedimentation-erosion
	$xtitle	$ytitle	12 0 1 3	 erosion-sedimentation rate
END




#SEDIMENT BUDGETS:
#volume in regions of interest:
#grdclip $tmp.sed.grd.tmp -G$tmp.sed.grd.tmp -Sb.01/NaN
if (-r $prj.INT) then 
	grd2xyz $tmp.sed.grd.tmp | outin $prj.INT | \
		awk -v dx=$dx -v dy=$dy '{if ($1=="1") {seds_in+=$4*dx*dy;}} END{printf("\nSediments in region 1: %.1f km3", seds_in);}'
endif
if (-r $prj.INT2) then 
	grd2xyz $tmp.sed.grd.tmp | outin $prj.INT2 | \
		awk -v dx=$dx -v dy=$dy '{if ($1=="1") {seds_in+=$4*dx*dy;}} END{printf("\nSediments in region 2: %.1f km3", seds_in);}'
endif
if (-r $prj.INT3) then 
	grd2xyz $tmp.sed.grd.tmp | outin $prj.INT3 | \
		awk -v dx=$dx -v dy=$dy '{if ($1=="1") {seds_in+=$4*dx*dy;}} END{printf("\nSediments in region 3: %.1f km3", seds_in);}'
endif
if (-r $prj.INT4) then 
	grd2xyz $tmp.sed.grd.tmp | outin $prj.INT4 | \
		awk -v dx=$dx -v dy=$dy '{if ($1=="1") {seds_in+=$4*dx*dy;}} END{printf("\nSediments in region 4: %.1f km3", seds_in);}'
endif
if (-r $prj.INT_ROSE) then 
	cat $prj.st | outin $prj.INT_ROSE | \
		awk -v dx=$dx -v dy=$dy '{if ($1=="1") {eros_in+=$4;topo_in+=$5;n_in++}} END{printf("\nMean Erosion in region INT_ROSE: %.1f m; mean topo: %.1f m", eros_in/n_in, topo_in/n_in);}'
endif




#2D CROSS SECTION:
set xscale = 0
set width_CS = `echo  2 \* $width + $separation  | bc -l`
if (-r $prj.pfl) then
echo plotting cross section...
set height_CS = 3.5
set vert_shift_CS = `echo -$height_CS \- $separation | bc -l`
set horz_shift_CS = `echo  -$width - $separation  | bc -l`
set dmax   	= `awk '{dmax=$3;} END{print dmax}' $prj.pfl`
psbasemap -JX$width_CS/$height_CS -R0/$dmax/$zmin/$zmax \
	-B$ticklabels\:"distance (km)":/$zticklabels\:"z (km)":SeW \
	-O -K -X$horz_shift_CS -Y$vert_shift_CS >> $ps
source $tisc_dir/script/tisc.common.cross_section.gmt.job
set xscale = `echo $width_CS / 2 | bc -l`
endif

endif

#PALETAS DE COLORES:
psscale -C$tmp.topo_palet.tmp 	-D$xscale/-1.5/$width_CS/.3h -B:"Elevation (m)": -L -O -K >> $ps
#psscale -C$tmp.eet_cpt.tmp 	-D$xscale/-3.1/$width_CS/.3h -B/:"EET": -O -L -K	>> $ps
#psscale -C$tmp.eros_cumul_cpt.tmp 	-D$xscale/-3.0/$width_CS/.3h -B:"Cumul. Sedimentation-Erosion (km)": -L -O -K >> $ps
psscale -C$tmp.erosrate_cpt.tmp -D$xscale/-3.0/$width_CS/.3h -B:"Erosion-Sedimentation rate (mm/yr)": -L -O -K >> $ps
#psscale -C$tmp.subs_cpt.tmp 	-D$xscale/-6.5/$width_CS/.3h -B/:"Subsidence (m)": -L -O -K >> $ps
psscale -C$tmp.rain_cpt.tmp 	-D$xscale/-4.5/$width_CS/.3h -B:"Precipitation (m/yr)": -L -O -K >> $ps

pstext	-JX -R -O -G0 -K -N -Y-1.2 <<END >> $ps 
	$xf	$y0	9 0 1 3	TISC software
END

#ROSE DIAGRAM:
if (-r $prj.INT_ROSE) then 
	cat $prj.xyw | outin $prj.INT_ROSE | awk '($6!="L" && !($10==$4 && $9==$3) && $1==1 && $4>.2){print $4, 90-atan2($9-$3,$8-$2)*180/3.1415927}' | \
		psrose -S1.5n -R0/1/0/360 -A45 -D -X0 -Y18.5 -Ba.5g.5/a45g45NE:."river orientation": -O -K -G220 >> $ps
endif


pstext -JX -R -O <<END>> $ps
END
endif



rm -f  $tmp.*.tmp

